Working with timestamps and bursts

This notebook is part of a tutorial series for the FRETBursts burst analysis software.

In this notebook we show how to access different streams of timestamps, burst data (i.e. start and stop times and indexes, counts, etc...). These operations are useful for users wanting to access or process bursts data and timestamps in custom ways. For a complete tutorial on burst analysis see FRETBursts - us-ALEX smFRET burst analysis.

Load data

We start by loading the data, computing background and performing a standard burst search:



In [ ]:

    
from fretbursts import *
sns = init_notebook()



In [ ]:

    
filename = "./data/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5"
d = loader.photon_hdf5(filename)
loader.alex_apply_period(d)
d.calc_bg(bg.exp_fit, time_s=30, tail_min_us='auto', F_bg=1.7)
d.burst_search()

Getting the timestamps

To get the timestamps arrays for a given photon stream we use Data.get_ph_times. Here a few example:



In [ ]:

    
ph = d.get_ph_times()                              # all the recorded photons
ph_dd = d.get_ph_times(ph_sel=Ph_sel(Dex='Dem'))   # donor excitation, donor emission
ph_d = d.get_ph_times(ph_sel=Ph_sel(Dex='DAem'))   # donor excitation, donor+acceptor emission
ph_aa = d.get_ph_times(ph_sel=Ph_sel(Aex='Aem'))   # acceptor excitation, acceptor emission

This are streams of all timestamps (both inside and outside the bursts). Similarly, we can get "masks" of photons for each photon stream using Data.get_ph_mask:



In [ ]:

    
mask_dd = d.get_ph_mask(ph_sel=Ph_sel(Dex='Dem'))   # donor excitation, donor emission
mask_d = d.get_ph_mask(ph_sel=Ph_sel(Dex='DAem'))   # donor excitation, donor+acceptor emission
mask_aa = d.get_ph_mask(ph_sel=Ph_sel(Aex='Aem'))   # acceptor excitation, acceptor emission

Masks are arrays of booleans (True or False values) which are True when the corresponding photon is in the stream. Note that all masks have same number of elements as the all-photons timestamps array:



In [ ]:

    
ph.size, mask_dd.size, mask_d.size, mask_aa.size

Masks can be used to count photons in one stream:



In [ ]:

    
mask_d.sum()

and to obtain the timestamps for one stream:



In [ ]:

    
ph[mask_d]

Note that the arrays ph[mask_d] and ph_d are equal. This is an important point to understand.

Burst data

There are several fields containing burst data:

Start-stop:

Data.mburst: start-stop information for each burst

Counts:

Data.nd: donor detector counts during donor excitation
Data.na: acceptor detector counts during donor excitation
Data.naa: acceptor detector counts during acceptor excitation (ALEX only)
Data.nda: donor detector counts during acceptor excitation

FRET:

Data.E: Proximity Ratio (when uncorrected) or FRET efficiency (when corrected)
Data.S: "Stoichiometry" (ALEX only)

All previous fields are lists containing one element per excitation spot. In single-spot data, these lists have only one element which is accessed using the [0] notation:



In [ ]:

    
bursts = d.mburst[0]
nd = d.nd[0]
na = d.na[0]
naa = d.naa[0]
E = d.E[0]
S = d.S[0]

All previous variables are numpy arrays, except for bursts which is a Bursts object (see next section).

Burst start and stop

The start-stop burst data is in bursts (a variable of type Bursts, plural):



In [ ]:

    
bursts

Indexing bursts we can access a single burst:



In [ ]:

    
firstburst = bursts[0]
firstburst

The first two "columns" (both in bursts or firstburst) are the index of first and last timestamps (relative to the all-photons timestamps). The last two columns (start and stop) are the actual times of burst start and stop. To access any of these fields we type:



In [ ]:

    
bursts.istart



In [ ]:

    
firstburst.istart

Note that ph[firstburst.istart] is equal to firstburst.start:



In [ ]:

    
ph[firstburst.istart], firstburst.start

The same holds for stop:



In [ ]:

    
ph[firstburst.istop], firstburst.stop

Note that bursts is a Bursts object (plural, a bursts-set) and firstburst is a Burst object (singular, only one burst). You can find more info on these objects in the documentation:

Low-level burst search functions

Burst photon-counts

The variables nd, na, naa contains the number of photon in different photon streams. By default these values are background corrected and, if the correction coefficients are specified, are also corrected for leakage, direct excitation and gamma factor.

To get the raw counts before correction we can redo the burst search as follow:



In [ ]:

    
d.burst_search(computefret=False)
d.calc_fret(count_ph=True, corrections=False)

Note that if you select bursts, you also need to use computefret=False to avoid recomputing E and S (which by default applies the corrections):



In [ ]:

    
ds = d.select_bursts(select_bursts.size, th1=30, computefret=False)



In [ ]:

    
nd = ds.nd[0]      # Donor-detector counts during donor excitation
na = ds.na[0]      # Acceptor-detector counts during donor excitation
naa = ds.naa[0]    # Acceptor-detector counts during acceptor excitation
E = ds.E[0]        # FRET efficiency or Proximity Ratio
S = ds.S[0]        # Stoichiometry, as defined in µs-ALEX experiments



In [ ]:

    
nd

Note that the burst counts are integer values, confirming that the background correction was not applied.

Slice bursts in time bins

Here we slice each burst in fixed time bins.



In [ ]:

    
from fretbursts.phtools.burstsearch import Burst, Bursts



In [ ]:

    
times = d.ph_times_m[0]  # timestamps array

We start fusing bursts with separation <= 0 milliseconds, to avoid having overlapping bursts:



In [ ]:

    
ds_fused = ds.fuse_bursts(ms=0)
bursts = ds_fused.mburst[0]
print('\nNumber of bursts:', bursts.num_bursts)

Now we can slice each burst using a constant time bin:



In [ ]:

    
time_bin = 0.5e-3  # 0.5 ms



In [ ]:

    
time_bin_clk = time_bin / ds.clk_p

sub_bursts_list = []
for burst in bursts:
    # Compute binning of current bursts
    bins = np.arange(burst.start, burst.stop + time_bin_clk, time_bin_clk)
    counts, _ = np.histogram(times[burst.istart:burst.istop+1], bins)
    
    # From `counts` in each bin, find start-stop times and indexes (sub-burst).
    # Note that start and stop are the min and max timestamps in the bin,
    # therefore they are not on the bin edges. Also the burst width is not
    # exactly equal to the bin width.
    sub_bursts_l = []
    sub_start = burst.start
    sub_istart = burst.istart
    for count in counts:
        # Let's skip bins with 0 photons
        if count == 0:
            continue
            
        sub_istop = sub_istart + count - 1
        sub_bursts_l.append(Burst(istart=sub_istart, istop=sub_istop,
                                  start=sub_start, stop=times[sub_istop]))
        
        sub_istart += count 
        sub_start = times[sub_istart]
    
    sub_bursts = Bursts.from_list(sub_bursts_l)
    assert sub_bursts.num_bursts > 0
    assert sub_bursts.width.max() < time_bin_clk
    sub_bursts_list.append(sub_bursts)

The list sub_bursts_list has one set of sub-bursts per each original burst:



In [ ]:

    
len(sub_bursts_list)



In [ ]:

    
ds_fused.num_bursts

Each set of sub-bursts is a usual Bursts object:



In [ ]:

    
print('Sub-bursts from burst 0:')
sub_bursts_list[0]



In [ ]:

    
iburst = 10
print('Sub-bursts from burst %d:' % iburst)
sub_bursts_list[iburst]

Photon counts in custom bursts

When performing a burst search, FRETBursts automatically computes donor and acceptor counts (in both excitation periods). These quantities are available as Data attributes: nd, na, naa and nda (as described in Burst data).

When a custom bursts-set is created, like in the previous section in which we sliced bursts in sub-bursts, we may want to computed the photon counts in the various photon streams. Let consider, as an example, the following Bursts object:



In [ ]:

    
bursts = sub_bursts_list[0]
bursts

How do we count the donor and acceptor photons in these bursts?

First we need to prepare the masks for the various photon streams (as explained before):



In [ ]:

    
mask_dd = d.get_ph_mask(ph_sel=Ph_sel(Dex='Dem'))   # donor excitation, donor emission
mask_ad = d.get_ph_mask(ph_sel=Ph_sel(Dex='Aem'))   # donor excitation, acceptor emission
mask_aa = d.get_ph_mask(ph_sel=Ph_sel(Aex='Aem'))   # acceptor excitation, acceptor emission

Next, we use the function counts_ph_in_bursts:



In [ ]:

    
from fretbursts.phtools.burstsearch import count_ph_in_bursts



In [ ]:

    
counts_dd = count_ph_in_bursts(bursts, mask_dd)
counts_dd

counts_dd contains the raw counts in each burst (in bursts) in the Donor-emission during Donor-acceptor stream. Similarly, we can compute counts for the other photon streams:



In [ ]:

    
counts_ad = count_ph_in_bursts(bursts, mask_ad)
counts_aa = count_ph_in_bursts(bursts, mask_aa)

With these values, we can compute, for example, the uncorrected Proximity Ratio (PR):



In [ ]:

    
counts_ad / (counts_dd + counts_ad)

Note that the first value is nan (not-a-number) because there is a zero in the numerator, i.e. the first bursts has 0 donor counts during donor excitation.

Finally, remember that no correction (background, leakage, etc.) has been applied so far.